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X) . Abstract 

We use stochastic computer simulations to study the transport of a spherical cargo particle along 
a microtubule-like track on a planar substrate by several kinesin-like processive motors. Our newly 
developed adhesive motor dynamics algorithm combines the numerical integration of a Langevin 
equation for the motion of a sphere with kinetic rules for the molecular motors. The Langevin 
I part includes diffusive motion, the action of the pulling motors, and hydrodynamic interactions 

between sphere and wall. The kinetic rules for the motors include binding to and unbinding 
from the filament as well as active motor steps. We find that the simulated mean transport 
length increases exponentially with the number of bound motors, in good agreement with earlier 
results. The number of motors in binding range to the motor track fiuctuates in time with a 
Poissonian distribution, both for springs and cables being used as models for the linker mechanics. 
Cooperativity in the sense of equal load sharing only occurs for high values for viscosity and 
attachment time. 

PACS numbers: 82.39.-k,87.10.Mn,87.15.hj 
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I. INTRODUCTION 



Molecular motors play a key role for the generation of movement and force in cellular 
systems [1]. In general there are two fundamentally different classes of molecular motors. 
Non-processive motors like myosin 11 motors in skeletal muscle bind to their tracks only 
for relatively short times. In order to generate movement and force, they therefore have to 
operate in sufficiently large numbers. Processive motors like kinesin remain attached to their 
tracks for a relatively long time and therefore are able to transport cargo over reasonable 
distances. Indeed processive cytoskeletal motors predominantly act as transport engines for 
cargo particles, including vesicles, small organelles, nuclei or viruses. For example, kinesin-1 
motors make an average of 100 steps of size 8 nm along a microtubule before detaching from 
the microtubule [2, 3], and therefore reach typical run lengths of micrometers. 

However, for intracellular transport even processive motors tend to function in ensembles 
of several motors, with typical motor numbers in the range of 1-10 [4]. The cooperation 
of several motors is required, for example, when processes like extrusion of lipid tethers 
require a certain level of force that exceeds the force generated by a single motor [5, 6]. The 
cooperative action of several processive motors is also required to achieve sufficiently long 
run length for cargo transport [7], as transport distances within cells are typically of the 
order of the cell size, larger than the micron single motor run length. In this context the most 
prominent example is axonal transport, as axons can extend over many centimeters [8, 9]. 
Another level of complexity of transport within cells is obtained by the simultaneous presence 
of different motor species on the same cargo, which can lead to bidirectional movements and 
switching between different types of tracks [10, 11], and by exchange of components of the 
motor complex with the cytoplasm. 

Cargo transport by molecular motors can be reconstituted in vitro using so-called bead 
assays in which motor molecules are firmly attached to spherical beads that flow in aqueous 
solution in a chamber. On the bottom wall of the chamber microtubules are mounted along 
which the beads can be transported [1, 12, 13, 14]. This assay has been used extensively to 
study transport by a single motor over the last decade [12], but recently several groups have 
adapted it for the quantitative characterization of transport by several motors [14, 15]. If 
several motors on the cargo can bind to the microtubule, then the transport process continues 
until all motors simultaneously unbind from the microtubule. Based on a theoretical model 



2 



for cooperative transport by several processive motors, it was recently predicted that the 
mean transport distance increases essentially exponentially with the number of available 
motors [7]. Indeed these predictions are in good agreement with experimental data [14, 16]. 
However, both the theoretical approach and the experiments do not allow us to investigate 
the details of this transport process. A major limitation of the bead assays for transport by 
several motors is that the number of motors per bead varies from bead to bead and that only 
the average number of motors per bead is known [14, 15]. In addition, even if the number 
of motors on the bead was known, the number of motors in binding range would still be a 
fluctuating quantity. Recently two kinesin motors have been elastically coupled by a DNA 
scaffold and the resulting transport has been analyzed in quantitative detail [16]. However, 
it is experimentally very challenging to extend this approach to higher numbers of motors. 

One key property of transport by molecular motors is the load force dependence of the 
transport velocity. For transport by single kinesins, the velocity decreases approximately 
linearly with increasing load and stalls at a load of about 6 pN [3]. Thus, when the cargo 
has to be transported against a large force, the speed of a single motor is slowed down. 
However, if several motors simultaneously pull the cargo, they could share the total load. 
This cooperativity lets them pull the cargo faster. Assuming equal load sharing, one can show 
that in the limit of large viscous load force the cargo velocity is expected to be proportional 
to the number of pulling motors [7]. Indeed, this is one explanation that was proposed 
by Hill et al. to give plausibility to their results from in vivo experiments, which showed 
that motor-pulled vesicles move at speeds of integer multiples of a certain velocity [9]. In 
general, however, one expects that the total load is not equally shared by the set of pulling 
motors. The force experienced by each motor will depend on its relative position along the 
track and can be expected to fluctuate due to the stochasticity of the motor steps [17]. In 
addition, for a spherical cargo particle curvature effects are expected to play a role of the 
way force is transmitted to the different motors. Because of its small size, the cargo particle 
is perpetually subject to thermal fluctuations. This diffusive particle motion is also expected 
to affect the load distribution and depends on the exact height of the sphere above the wall 
due to the hydrodynamic interactions. 

In order to investigate these effects, here we introduce an algorithm that allows us to 
simulate the transport of a spherical particle by kinesin-like motors along a straight filament 
that is mounted to a plain wall. Binding and unbinding of the motor to the filament 
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can be described in the same theoretical framework as the reaction dynamics of receptor- 
hgand bonds [18, 19, 20]. Similarly as receptors bind very specifically to certain ligands, 
conventional kinesin binds only to certain sites on the microtubule. Thus from the theoretical 
point of view a spherical particle covered with motors binding to tracks on the substrate is 
equivalent to a receptor-covered cell binding to a ligand-covered substrate. This situation is 
reminescent of rolling adhesion, the phenomenon that in the vasculature different cell types 
(mainly white blood cells, but also cancer cells, stem cells or malaria-infected red blood 
cells) bind to the vessel walls under transport conditions [21]. Different approaches have 
been developed to understand the combination of transport and receptor-ligand kinetics in 
rolling adhesion. Among these Hammer and coworkers developed an algorithm that combines 
hydrodynamic interactions with reaction kinetics for receptor-ligand bonds [22]. Recently, 
we introduced a new version of this algorithm that also includes diffusive motion of the 
spherical particle [23]. Here, we further extend our algorithm to include the active stepping 
of motors {adhesive motor dynamics) . Simulation experiments with this algorithm provide 
access to experimentally hidden observables like the number of actually pulling motors, 
the relative position of the motors to each other, and load distributions. The influence of 
thermal fluctuations on the motion of the cargo particle is also influenced by the properties 
of the molecules that link the cargo to the microtubule. In our simulations various polymer 
models can be implemented and their influence can be tested directly. In general our method 
makes it possible to probe the effects of various microscopic models for motor mechanics on 
macroscopic observables that are directly accessible to experiments. 

The organization of the article is as follows. In the first part. Sec. II, we explain our model 
in detail. This is based on a Langevin equation that allows us to calculate the position and 
orientation of a spherical cargo particle as a function of time. In addition, we include rules 
that model the reaction kinetics of the molecular motors being attached to the cargo and 
comment on the different kinds of friction involved. We then explain how theoretical results 
for the dependence of the mean run length of a cargo particle on the number of available 
motors previously obtained in the framework of a master equations can be compared to the 
situation where only the total number of motors attached to the cargo sphere is known. 
We also briefly comment on the implementation of our simulations. In the results part. 
Sec. Ill, we first measure the mean run length and the mean number of pulling motors at low 
viscous drag and find good agreement with earlier results. We then present measurements 
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of quantities that are not accessible in earlier approaches, including the dynamics of the 
number of motors on the cargo that are in binding range to the microtubule. Finally, 
we consider cargo transport in the high viscosity regime and investigate how the load is 
distributed among the pulling motors. We find that cooperativity by load sharing strongly 
depends on appropriate life times of bound motors. In the closing part. Sec. IV, we discuss 
to what extend our simulations connect theoretical modeling with experimental findings. 
Furthermore, we give an outlook on further possible applications of the adhesive motor 
dynamics algorithm introduced here. 

II. MODEL AND COMPUTATIONAL METHODS 
A. Bead dynamics 

In experiments using bead assays for studying the collective transport behavior of kinesin 
motors one notes the presence of three very different length scales, namely the chamber 
dimension, the bead size and the molecular dimensions of kinesin and microtubule, respec- 
tively. The chamber dimension is macroscopic. The typical radius R of beads in assay 
chambers is in the micrometer range [13]. The kinesin molecules with which they are cov- 
ered have a resting length /q of about 80 nm [24] . Kinesins walk along microtubules which are 
long hollow cylindrical filaments made from 13 parallel protofilaments and with a diameter 
of about Hmt = 24 nm [25]. Thus the chamber dimensions are large compared to the bead 
radius, which in turn is large compared to the motors and their tracks. This separation of 
length scales allows us to model the microtubule as a line of binding sites covering the wall 
and means that the dominant hydrodynamic interaction is the one between the spherical 
cargo and the wall. For sufficiently small motor density, this separation of length scales also 
implies that we have to consider only one lane of binding sites. In the following we therefore 
consider a rigid sphere of radius R moving above a planar wall with an embedded line of 
binding sites as a simple model system for the cargo transport by molecular motors along a 
filament. 

For small objects like microspheres typical values of the Reynolds number are much 
smaller than one and inertia can be neglected (overdamped regime). Therefore, the hydro- 
dynamic interaction between the sphere and the wall is described by the Stokes equation. 
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Throughout this paper we consider vanishing external flow around the bead. Directional 
motion of the sphere arises from the pulling forces exerted by the motors. In addition the 
bead is subject to thermal fluctuations that are ubiquitous for microscopic objects. Tra- 
jectories of the bead are therefore described by an appropriate Langevin equation. For the 
sake of a concise notation we introduce the six- dimensional state vector X which includes 
both the three translational and the three rotational degrees of freedom of the sphere. The 
translational degrees of freedom of X refer to the center of mass of the sphere with respect 
to some reference frame (cf. Fig. 1). The rotational part of X denotes the angles by which 
the coordinate system fixed to the sphere is rotated relative to the reference frame [26]. 
Similarly, F denotes a combined six- dimensional force/torque vector. 

With this notation at hand the appropriate Langevin equation reads [27, 28] 



Here, M is the position-dependent 6x6 mobility matrix. As we consider no-slip boundary 
conditions at the wall, M depends on the height of the sphere above the wall in such a way 
that the mobility is zero when the sphere touches the wall [29]. Thus, the hydrodynamic 
interaction between the sphere and the wall is completely included in the configuration 
dependence of the mobility matrix. The last term in Eq. (1) is a Gaussian white noise term 
with 



with Boltzmann's constant ks- The second equation represents the fluctuation-dissipation 
theorem illustrating that the noise is multiplicative due to the position-dependance of M. 
Thus we also have to define in which sense the noise in Eq. (1) shall be interpreted. As 
usual for physical processes modeled in the limit of vanishing correlation time we choose 
the Stratonovich interpretation [30]. However, Eq. (1) is written in the Ito version marked 
by the super-index / for the noise term. The Ito version provides a suitable base for the 
numerical integration of the Langevin equation using a simple Euler scheme. The gradient 
term in Eq. (1) is the combined result of using the Ito version of the noise and a term that 
compensates a spurious drift term arising from the no-slip boundary conditions [26, 27]. 



X = MF + ksTS/M + gl 



(1) 



(g[) = 0, {glgl,) = 2kBTM5{t-t'), 



(2) 
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FIG. 1: A single sphere of radius R and surface separation h from a planar wall. The translational 
coordinates R of the sphere are given relative to a reference frame that is fixed to the wall. The 
sphere is pulled by one molecular motor that is attached to the surface at position r measured with 
respect to the center of the sphere. The bead is subject to the motor force with x-component 
Fm,x- In addition, an external force Ft acts parallel to the filament, typically arising from an 
optical trap. The force unbalance between Fm^x and Ft leads to the bead velocity U. The motor 
with resting length is firmly attached to the bead and can bind to and unbind from a MT and 
moves with velocity Vm- X denotes the angle between motor and MT. 

For our numerical simulations we discretize Eq. (1) in time and use an Euler algorithm 
which is of first order in the time step At 

A'Xt = FAt + kBTVMAt + g{At) + 0{Af) (3) 

where g{At) has the same statistical properties as gf from above. In order to compute the 
position-dependent mobility matrix M we use a scheme presented by Jones et al. [29, 31] 
that provides accurate results for all values of the height z. A detailed description of the 
complete algorithm including the translational and rotational update of the sphere can be 
found in Ref. [26]. 

B. Motor dynamics 

In our model the spherical cargo particle is uniformly covered with Ntot motor proteins. 
These molecules are firmly attached to the sphere at their foot domains. Consequently, 
Ntot is a constant in time. The opposite ends of these molecules (their head domains) can 
reversibly attach to the microtubule (MT), which is modeled as a line of equally spaced 
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binding sites for the motors covering the wall. A motor that is bound to the MT exerts a 
force and a torque to the cargo particle. If Vh and rj are the positions of the head and foot 
domains of the motor, respectively, then the force by the motor is given by: 

where the absolute value of the motor force is given by the force extension relation 
F{x) for the motor protein. Throughout this article we consider two variants of the force 
extension curve for the polymeric tail of the motor. The harmonic spring model reads 

F(x) = k(x-/o). (5) 

This means, a force is needed for both compression and extension of the motor protein. 
Actually, it was found that kinesin exhibits a non-linear force extension relation [32], with 
the spring constant varying between k = 0.2 • 10~^ N/m for small extensions and k = 
0.6 ■ 10~^ N/m for larger extensions [32, 33]. For extensions close to the contour length 
the molecule becomes infinitely stiff [34]. For small extensions, however, the harmonic 
approximation works well. Alternatively, we consider the cable model 

F{x) = k{x - loMx - /o), e(a:) = { ' ^ V (6) 

else 

In the cable model force is only built up when the motor is extended above its resting length 
Iq. In the compression mode, i. e. when the actual motor length is less then the resting 
length, no force exists. The cable model can be seen as the simplest model for a flexible 
polymer. 

Besides the force each motor attached to the MT also exerts a torque on the cargo particle. 
With f being the position of the motor foot relative to the center of the sphere (cf. Fig. 1), 
this torque reads 

= r X F„. (7) 

The combined force/torque vector F = (F^, T^)"^ enters the Langevin algorithm Eq. (3). 
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In addition to the firm attachment of the motors to the sphere, each motor can in principle 
reversibly bind and unbind to the MT. We model these processes as simple Poissonian rate 
processes in similar a fashion as it is done for modeling of formation and rupture of receptor 
ligand bonds in cell adhesion (e. g., [22]). The motor head is allowed to rotate freely about its 
point of fixation on the cargo. The head is therefore located on a spherical shell with radius 
given by the motor resting length Iq. However, in contrast to the anchorage point on the cargo 
particle the head position of the motor is not explicitely resolved by the algorithm. In order 
to model the binding process we introduce a capture length tq. A motor head can then bind 
to the MT with binding rate iTad whenever the spherical shell of radius Iq and thickness ro 
around the motor's anchorage point has some overlap with a non-occupied binding site on the 
MT. The MT's binding sites are identified by the tubulin building block of length 6 = 8 nm, 
so we choose tq = 6/2. Note that binding rate and binding range are complementary 
quantities and that a more detailed modeling of the binding process would have to yield 
appropriate values for both quantities. If the overlap criterion is fulfilled within a time 
interval At, the probability for binding pon within this time step is pon = 1 — exp(— TTaj^At). 
With a standard Monte-Carlo technique it is then decided whether binding occurs or not: 
a random number rand is drawn from the uniform distribution on the interval (0, 1) and in 
the case Pon > rand binding occurs. If the overlap criterion is fulfilled for several binding 
sites, then using the Monte-Carlo technique it is first decided whether binding occurs and 
then one of the possible binding sites is randomly chosen. 

Each motor bound to the MT can unbind with escape rate e. In single motor experiments 
it was found that the escape rate increases with increasing motor force [35]. This force 
dependence can be described by the Bell equation [7, 38] 

e = eoexp{EjFa), (8) 

with the unstressed escape rate eo and the detachment force F^. Because the details of forced 
motor unbinding from a filament are not known, here we make the simple assumption that 
the unbinding pathway is oriented in the direction of the tether. We set e = eo whenever 
the motor is compressed. 

The major conceptual difference between a motor connecting a sphere with a MT and 
a receptor-ligand bond is that a motor can actively step forward from one binding site to 
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Parameter 


typical value 


meaning 


reference 


R 


1 fim 


bead radius 




eo 


1 


unstressed escape rate 


[7] 




5 s"^ 


binding rate 


[6,7] 


Fd 


3 pN 


detachment force 


[7, 35] 


K 10"^ 


, ..10-^...10"3 


N/m motor spring constant 


[33, 34, 36] 


5 


8 nm 


kinesin step length 


[2] 


vo 


1 /im/s 


maximum motor velocity 


[7] 


A° := vo/6 


125 


forward step rate 




ro 


5/2 


capture radius 




Fs 


5 ... 6 ... 8 pN 


stall force 


[3] 


lo 


50,65,80 nm 


(resting) length 


[24, 37] 




24 nm 


microtubule diameter 


[25] 



TABLE I: Parameters used for adhesive motor dynamics. For ambient temperature we use T = 
293 K, for viscosity rj = 1 mPas (if not otherwise stated). If a range is given, then figure in bold 
face denotes the value used in the numerical simulations. 



the next with step length 6 (the length of the tubulin unit). The mean velocity Vq of an 
unloaded kinesin motor is about vq = 1 /im/s depending on ATP concentration [3]. If the 
motor protein is mechanically loaded with force opposing the walking direction, the motor 
velocity Vm is decreased. For a single kinesin molecule that is attached to a bead on which a 
trap force Rt pulls, the velocity was found experimentally to decrease approximately linearly 
[3, 39]: 

Vm =Vo(^l-^^, 0<Rt< Rs, (9) 

with the stall force Rs and the trap force Rt acting antiparallel to the walking direction. 
Different experiments have reported stall forces between 5 and 8 pN. Changes in this range 
are not essential for our results and therefore we use the intermediate value Rb = Q pN. 
If the force is higher than the stall force, kinesin motors walk backwards with a very low 
velocity [40], which we will neglect in the following. Finally, for assisting forces, i.e. if the 
motor is pulled forward, the effect of force is relatively small [40, 41]. In order to derive 
an expression similar to Eq. (9) for our model, we have to identify the proper term that 
replaces Rt in Eq. (9). First, we rewrite Eq. (9) {Rs — Rt), with some internal 

motor mobility coefficient /i^ := Vq/Rs- This version of Eq. (9) allows us to interpret the 
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m,x I 



FIG. 2: Force velocity relation for a single motor: velocity Vm as a function of load |-Fm,x| according 
to Eq. (10) with maximum velocity vq and stall force K,. 



motor head as an over-damped (Stokesian) particle that constantly pulls with the force Fg 
against some external load Ft resulting in the effective velocity Vm- According to Eq. (4) a 
motor pulls with force on the bead, so we can identify the "load" to be —F^^r^, where 
is the x-component of and the minus sign accounts for Newton's third law {actio = 
reactio). Thus, we obtain the following piecewise linear force velocity relation for the single 
motor (see also Fig. 2): 



Vm = Vo< 



I -^m , X I 







if 



< 



if < F„ ■ < 



(10) 



if 



F ■ p > F 



where e^; is the walking direction of the motor, see Fig. 1. Thus, if the motor pulls antiparallel 
to its walking direction on the bead, it walks with its maximum speed vq. If it is loaded with 
force exceeding the stall force Fg, it stops. For intermediate loadings the velocity decreases 
linearly with load force. Eq. (10) defines the mean velocity of a motor in the presence of 
loading force. We note that our algorithm also allows us to implement more complicated 
force- velocity relations and is not restricted to the piecewise linear force- velocity relation. 
It is used here because we do not focus on a specific kinesin motor and because it is easy 
to implement in the computer simulations. In practice the motor walks with discrete steps 
of length S. In the algorithm we account for this by defining a step rate := Vm/^- The 
decision for a step during a time interval At is then made with the same Monte-Carlo 
technique used to model the binding and unbinding process of the motor head. A step is 
rejected if the next binding site is already occupied by another motor (mutual exclusion) 
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[42]. 



C. Bead versus motor friction 

In principle, the velocity of a single motor Vm pulling the sphere and the component in 
walking direction of the sphere's velocity U (see Fig. 1) are not the same. For an external 
force Ft (against walking direction) in the pN range acting on the sphere and -Fm,x being 
the motor force in walking direction, the bead velocity U is given hj U = fi^J:^{Fm,x — Ft),. 
Here, /xjf^ denotes the corresponding component of the mobility matrix M of the sphere 
(cf. Eq. (1)) evaluated at the height of the sphere's center with super indices tt referring 
to the translational sector of the matrix and sub indices xx referring to the responses in 
x-direction. On the other hand from Eq. (10) it follows that the motor head moves with 
velocity Vm = fim.{Fs — Fm,x)- Only in the stationary state of motion the two speeds are 
equal, U = Vm, and we obtain the force with which the motor pulls (in walking direction) 
on the bead: 

Fn.,x = -^Ft + -^F,. (11) 
A'-m + ^^xx f^m I f^xx 

Thus, if the internal friction of the motor is large compared to the viscous friction of the 
sphere, i.e., l//im ^ ^^e second term in Eq. (11) can be neglected and one has 

Fm,x ~ Ft. That means, only the trap force pulls on the motor. If fi^ ~ f^xx both terms in 
Eq. (11) are of the same order of magnitude. Then, both external load Ft and the friction 
force on the bead will influence the motor velocity. Experimentally, these prediction can 
be checked in bead assays by varying the viscosities of the medium (e.g., by adding sugar 
like dextran or FicoU [43]). Numerically, we can vary rj in the adhesive motor dynamics 
algorithm. 

When several motors are simultaneously pulling, they can cooperate by sharing the load. 
Assuming the case that n motors are attached to the MT which equally share the total load, 
then the force in x-direction exerted on the bead is nFm.x{n) with Fm,x being again the 
individual motor force. In the stationary state with U = Vm we have (with external load 
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Ft = 0): 



Fn.An) = F,. (12) 

The mean bead velocity U{n) with n equally pulling motors is U{n) = nfi^J:^Fm,x{n). Thus, 
in general, the velocity of the beads will increase with increasing number n of pulling motors. 
However under typical experimental conditions in vitro, where bead movements are probed 
in aqueous solutions, the internal motor friction l/fim dominates over the viscous friction 
of the bead, l/f^c^, and the velocity becomes independent of the number of motors. Only 
if the bead friction becomes comparable to the internal motor friction, the velocity exhibits 
an appreciable dependence on the motor number. This dependence can be illustrated by 
considering the ratio of U{n) and U{1): 

Uil) n/ig^ + flm n/ fim + 1/ ' 

In the opposite limit, n/fim <^ V/^xx; i-^-, when the viscous bead friction is very large and 
dominates over the internal motor friction, Eq. (13), leads to U{n) ^ nU{l), and the bead 
velocity increases linearly with n [7]. 



D. Vertical forces 



We note that, although we are mainly interested in the x-component of the motor force 
Fm, i-e. the component parallel to the microtubule, which enters the force- velocity relation, 
the motor force also has a component perpendicular to the microtubule, see Fig. 1. This 
force component tends to pull the bead towards the microtubule and thus to the surface, 
whenever the bead is connected to the microtubules by a motor. This force is balanced by 
the microtubule repelling the bead. Additionally, if the bead touches the filament or the 
wall, diffusion can only move the bead away from the wall. In case of the full spring model, 
compressed motors can also contribute a repulsion of the bead from the wall. If viscous 
friction is strong, the normal component of the force arises mainly from the microtubule. In 
that case, the distance h between the bead and the microtubule is very small, (h) ^ 0. For 
small viscous friction, thermal fluctuations play a major role and lead to non-zero distances 
between the bead and the microtubule, as discussed further in Sec. IIIB. 
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All parameters used for the adhesive motor dynamics simulations together with typical 
values are summarized in Tab. I. For the numerical simulations we non-dimensionalize all 
quantities using R for the length scale, 1/eo for the time scale and the detachment force Fa 
as force scale. 

E. Master equation approach 

When no external load is applied a motor walks on average a time 1/eo before it detaches 
and the cargo particle might diffuse away from the MT. When several motors on the cargo 
can bind to the MT the mean run length dramatically increases as was previously shown 
with a master equation approach [7]. For the sake of later comparison to our simulations 
we briefly summarize some of these results in the following. 

Let Pi be the probability that i motors are simultaneously bound to the MT with 
i = 0, . . . , Nm and Nm being the maximum number of motors that can bind to the MT 
simultaneously. Assuming the system of Nm motors is in a stationary state and the total 
probability of having i = 0, . . . , N^, motors bound is conserved then the probability flux 
from one state to a neighboring state is zero. This means that the probability Pj of having 
i bound motors can be calculated by equating forward and reverse fluxes 

{Nm - i)TradPi = + l)eP^+u ^ = 0, . . . , iV^ - 1, (14) 

where it is assumed that the escape rate e is a constant with respect to time. The solution 
to Eq. (14) is given by 



P^= I " ^—^^ 2 = 0,. ..,A^^. (15) 



e + vr 



The probability that i motors are simultaneously pulling under the condition that at least 
one motor is pulling is -Pi/(1 — -Pq) for ^ = 1, • • • , A^m- Then, the mean number of bound 
motors Ni, (given that at least one motor is bound) is [7, Eq. [13]] 
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FIG. 3: Illustration of the area fraction a^, of the sphere (cut and placed besides the sphere) on 
which motor proteins can reach the MT (thin cylinder), a;, depends in a geometrical fashion on 
the minimum distance h between the sphere and the MT and on the resting length ^o- 

The effective unbinding rate eg//, i. e., the rate with which the system reaches the unbound 
state, is determined from ee//(l — Pq) = T^adPo- This quantity can also be identified with the 
inverse of the mean first passage time for reaching the unbound state, when starting with 
one motor bound [7]. If the medium viscosity is small, i.e., similar to that of water, and no 
external force is pulling on the bead, we assume that the velocity of the bead U does not 
depend on the number of pulling motors. The mean run length, that is the mean distance 
the cargo is transported by the motors in the case that initially one motor was bound, is 
then the product of mean velocity U and mean lifetime (l/eg//) [7, Eq. [14]]: 



For kinesin-like motors at small external load with iTad ^ e this expression can be ap- 
proximated by {Axb) ~ {U/Nm.e){Tiad/^)'^"^~^, i-e., the mean run length grows essentially 
exponentially with N^- In the stationary state the bead velocity U and the motor velocity Vm 
are equal. For no external load and small viscous friction on the bead one can approximate 
e eo in Eq. (16) and Eq. (17). 

F. Mean run length for a spherical cargo particle 

In contrast to the master equation model in which one fixes the maximal number of 
bound motors Nm, in the computer simulations only the total number Ntot of motors on the 
spherical cargo particle is fixed. A similar situation arises in experiments where only the total 




(17) 
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amount of molecular motors on the sphere is measured [14] (but not in experiments with 
defined multi-motor complexes [16]). If Ab is the area on the sphere's surface that includes 
all points being less than Iq apart from the MT (cf. Fig. 3), we expect on average Uf, = Ntotih 
motors to be close enough to the MT for binding, with the reduced area := Ai, / {4:TT R^) . 
While the master equation model assumes a fixed Nm-, in the simulations the maximum 
number of motors that can simultaneously bind to the MT is a fiuctuating quantity about 
the mean value n^. In the simulations the motors are uniformly distributed on the cargo, 
thus the probability distribution function P{k) for placing k motors inside the above defined 
area fraction is a binomial distribution. As we have Iq <^ i?, a?, is small and P{k) is well 
approximated by the Poissonian probability distribution function 
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P{k,nb) = ^exp(-nb). 



nb = Ntotdb- 



(18) 



Thus, given a fixed number Ntot of motors on the sphere the number of motors that are 
initially in binding range to the MT might be different from run to run. In addition, 
because of thermal fiuctuations and torques that the motors may exert on the cargo particle 
the relative orientation of the sphere changes during a simulation run. With the change of 
orientation also the number of motors in binding range to the MT is not a constant quantity 
during one run. 

In order to make simulation results for the mean run length and the mean number of 
bound motors comparable with the theoretical predictions Eq. (17) and Eq. (16), respec- 
tively, we have to average over different Nm- Neglecting fluctuations in the number of 
motors that are in binding range to the MT during one simulation run we perform the aver- 
age with respect to the Poisson distribution, Eq. (18). Averaging the mean walking distance 
{Axb){Nm) from Eq. (17) over all Nm with weighting factors given by Eq. (18) we obtain 
the following expression (rif, = abNtot)'- 
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(19) 



We note that during the initialization of each simulation run we place one motor on the 
lower apex of the sphere and then distribute the other Ntot — 1 motors uniformly over 
the whole sphere (see appendix B for a detailed description of the procedure). For this 
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reason P{Nm — l,nb) denotes the probability of having in total N„i motors inside the area 
fraction at,. Furthermore, we introduced a cutoff Nm.max of maximal possible motors {Ntot 
is obviously an upper limit for Nm,max)- In the hmit N^^rnax oo we have {{Axb))poisson = 

In a similar way we can calculate the Poisson-averaged mean number of bound motors 
{Nb)poisson- For this it is important to include the correct weighting factor [44]. From 
an ensemble of many simulation runs those with larger run length contribute more than 
those with smaller run length. For simulation runs the mean number of bound motors 
is obtained as {Nb)sim = Xlj ^i'^(^)/ Xli where ti is a period of time during which n{i) 
motors are bound and the sum is over all such periods of time. Assuming the bead velocity 
f/ to be a constant, the time periods ti can also be replaced by the run lengths Axb,i during 
ti. Picking out all simulation runs with a fixed N^, their contribution to the sum is the 
mean number of bound motors A^;, times the total run length of beads with given A^^- The 
latter is the mean run length (Axb) Nm times the number of simulation runs with the given 
Nm (for sufficiently large A^). Clearly, the fraction of runs with given Nm is the probability 
P{Nm — l,nb) introduced in Eq. (18). Consequently, we obtain again with a truncation of 
the sum at some Nm,max < Ntot 

""x^^^ PfM ^^rf^r ^ Et=T PjNm - l,nb){Axb)N„.Nb{Nm) 

{Nb)po^sson= PiNm)NbiNm) = ^„ ■ (20) 

N^=l 2^N^=l P{Nm-l,nb){Axb}N^ 

Here we introduced the probability P{n) for having n motors in binding range to the MT 
when picking out some cargo particle from a large ensemble of spheres. Explicitely, this 
probability is given by 

n-l 

p. . _ £((1 + ^^'^/^r - 1) Ve-^ _ ((1 + n^Jeor - l)nVln\ .31) 
G. Computer simulations 

We use the following procedure for the computer simulations. In each simulation run the 
sphere is covered with Ntot motors. Initially, one motor, located at the lowest point of the 
sphere, is attached to the microtubule such that the distance of closest approach h between 
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the sphere and the microtubule is given by the resting length of the motor, i. e., h = Iq. The 
other [Ntot — 1) motors are uniformly distributed on the sphere's surface (cf. Sec. B). When 
the motor starts walking, it pulls the sphere closer to the MT because there is a ^-component 
in the force exerted on the sphere by the motor stalk (which is strained after the first step). 
Then, other motors can bind to the MT. The system needs some time to reach a stationary 
state of motion, so initially the motor velocity Vm and the bead velocity U are not the same 
(for reasons of comparison a fixed initial position is necessary; other initial positions have 
been tested but initialization effects were always visible). In principle a simulation run lasts 
until no motor is bound. For computational reasons, each run is stopped after 2-10^ s 
(which is rarely reached for the parameters under consideration). For each run quantities 
like the mean number of bound motors Nh, the walking distance Axb and the mean distance 
of closest approach (h) between sphere and MT are recorded. 

III. RESULTS 

A. Single motor simulations 

In Sec. II B we defined a force- velocity relation for the single motor. In this section we 
perform simulation runs with a single motor, i.e., Ntot = 1, and measure the effective force 
velocity relation by tracking the position of the sphere. Inserting Eq. (11) into Eq. (10) 
provides a prediction for the velocity of a bead subject to one pulling motor and an external 
trap force Ft. In Fig. 4a this prediction is shown for three different values of viscosity 
together with the actual measured velocities during the simulations. More precisely, we 
measured the mean velocity of the bead and the motor obtained from a large number of 
simulation runs (to avoid effects resulting from the initial conditions we first allowed the 
relative position/orientation of bead and motor to "equilibrate" before starting the actual 
measurement). The mean velocity is then given as the total (summed up over all simulation 
runs) run length divided by the total walking time. The good agreement between the 
numerical results and the theoretical predictions provides a favorable test to the algorithm. 
At ?7 = 1 mPas (the viscosity of water), friction of the bead has almost no influence on 
the walking speed. At hundred times larger viscosities, however, bead friction reduces the 
motor speed to almost half of its maximum value already at zero external load. Although 
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(a) (b) 

FIG. 4: (a) Measured force velocity relation of a single motor (with Iq = 80 nm) pulling a sphere 
of radius R = 1 fim for three different viscosities 77 = 1,10,100 mPas. Shown is the relation 
according to Eq. (9), the actual measured force velocity relation of the motor head and the bead 
center, respectively, and the theoretical prediction according to Eq. (10) and Eq. (11). (b) The 
measured force velocity relation for r] = 1 mPas is shown where in Eq. (10) not |-FTn,a:| but ||Fm|| is 
used. The dotted line emphasizes the linear decrease of the velocity. The negative velocity of the 
bead at large Ft results from thermal fluctuations. Fluctuations against walking direction increase 
the escape probability. In case of escape they cannot be compensated by fluctuations in walking 
direction. (Numerical parameters: At = 10~^, number of runs iV = 2 - 10^ - 9- 10^.) 



the velocities of the motor and the bead are expected to be equal, Fig. 4a shows that the 
motor is slightly faster than the bead. This is a result of the discrete steps of the motor 
and can be considered as a numerical artefact: at the moment the motor steps forward the 
motor stalk is slightly more stretched (loaded) than before the step, therefore, the escape 
probability is increased. The result of unbinding at the next time step would then be that the 
bead moved a distance 6 less than the motor. For loads close to the stall force the observed 
velocity is slightly larger than the prediction, which is a result of thermal fluctuations of the 
bead in combination with the stepwise linearity of the force velocity relation: a fluctuation 
in walking direction slightly reduces the load on the motor, thus increasing the step rate, 
whereas fluctuations against walking direction lead to zero step rate. 

It was observed by Block et al. that vertical forces on the bead (i. e., in 2;-direction) also 
reduce the velocity of the motor [41]. But the same force that leads to stall when applied 
antiparallel to the walking direction has a rather weak effect on the motor velocity when 
applied in ^-direction. Using the absolute value of the total force of the motor ||Fm|| in 
Eq. (10) instead of its x-component |-Fm,2:|, we measure a force velocity relation as shown in 
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FIG. 5: Mean run length (Ax^) of a bead pulled by a single motor as a function of an external force 
on the bead Ft and for three different viscosities = 1, 10, 100 mPas. The lines give the theoretical 
predictions according to Eq. (22) assuming an angle of 60° between the motor and the MT. For 
comparison also the theoretically predicted (A3;{,)-curve for x = is shown (double dotted line). 



; 7/ = lOOmPas ' + - v = lOnmPas, x = 60° 

■ 77 = 10mPas ^-X--: >? = lOmPas, X = 60° ■ 
/; = ImPas I X I = imPas, \ = 60° ■ 
q = ImPas, X = 0° ■ 




Fig. 4b. Again, the velocity decreases essentially linearly with applied external force, but 
stalls already at around Ft ~ Fs/2 because of vertical contributions of the force ||F,„||. As 
the effect of vertical loading reported in Ref. [41] seems to be much weaker than that shown 
in Fig. 4b, we reject this choice of force velocity relation. 

From the simulations carried out for Fig. 4a we can also obtain the mean run length 
(Axf,) for a single motor as a function of external load. The results are shown in Fig. 5. 
Using Eq. (10) and the Bell equation, Eq. (8), we obtain 

^ e "6oexp(||F„||/F,)- ^''^ 

The numerical results shown in Fig. 5 fit very well to the theoretical prediction of Eq. (22) 
when assuming the angle x between the motor and the MT to be x = 60°. The angle x 
depends on the bead radius i?, the resting length /q [45] and the polymer characteristics of 
the motor protein, e. g., its stiffness k. 



B. Run length for several motors pulling 

We now measure the run length distributions and the mean run length as a function of 
motor coverage Ntot- For motors modeled as springs according to Eq. (5) with two different 
resting lengths /q = 50, 65 nm the run length distributions are shown in Fig. 6a,b. For each 
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FIG. 6: Distribution of run lengths Axf, in semi-logarithmic scale for different values of motor 
coverage Ntot- The motor protein is modeled as a harmonic spring according to Eq. (5). (a) 
Resting length of the motor protein Iq = 50 nm. (b) Resting length of the motor protein Iq = 65 nm. 
(Numerical parameters: time step At = 10~^/eo, number of simulation runs N ~ 10^.) 

value of Ntot the run length was measured about = 10^ times. The simulations turn 
out to be very costly, especially for large Ntot as the mean run length increases essentially 
exponentially with the number of pulling motors (cf. Eq. (17)). From Fig. 6 we see that the 
larger Ntot, the more probable large run lengths are, resulting in distribution functions that 
exhibit a flatter and flatter tail upon increasing Ntot- Fig. 6 nicely shows that the shape of 
the distributions depends not only on the total number of motors Ntot attached to the sphere 
but also on the resting lengths. Given the same Ntot "we can see that longer run lengths are 
more probable for the longer resting length Iq = 65 nm shown in Fig. 6b. This can simply 
be explained by the fact that the larger the motor proteins, the larger is the area fraction 
ttf, introduced in Fig. 3 and therefore the more motors are on average close enough to the 
microtubule to bind. 

Fig. 7a shows the mean run length as a function of Ntot as obtained by numerical simu- 
lations of the transport process (points with error bars). For the motor stalk three different 
values of the resting length Iq = 50, 65, 80 nm are chosen and both the full-spring and the 
cable model are applied for the force extension relation. Similarly to what we have already 
seen for the run length distributions in Fig. 6, the larger the resting length Iq the more mo- 
tors can simultaneously bind for given Ntot, and therefore the larger is the mean run length. 
Furthermore, Fig. 7a also demonstrates that it makes a clear difference whether the motor 
stalk behaves like a full harmonic spring or a cable. If the motor protein behaves like a cable 



21 



10- 



E 10'' r 



10' 



10' 



10" 



la = 50 nm, spring i — I — i 

l(, = 50 nm, cable ■— X--' 
k = 65 nm, spring ■■ 

k = 65 nm, cable :■■ ■■: ..■® 
la = 80 nm, spring < ■ < ,.®' 

(ii = 80 nm, cable i— ©— i g,-'' 

fit for area-fraction 

80 160 240 320 400 
total number of motors N,, 

fa) 



® 



X.- 



.X- 

...+■ 



/,) = 50 nm, spring 

If, = 80 nm, spring ^ X 

la = 50 nm, cable 

la = 80 nm, cable D - 



X'" T 

* 



- * 



X 



T * 



480 560 



80 160 240 320 400 480 560 
total number of motors N,„, 



(b) 



FIG. 7: (a) Mean run length {Axb) (data points with error bars) as a function of motors on the 
bead Ntot obtained from adhesive motor dynamics. The hues are fits of Eq. (19) with respect to 
the area fraction ah- (b) Mean number of bound motors Ni, (data points with error bars). The 
Unes are the values obtained from the Poisson-averaeed mean number of bound motors {^b) Poisson 
in Eq. (20) using for the fit value from (a). (Parameters: -Kad = 5, eq = 1, = 125, Ai = 10~^, 
N ~ 10"^.) 



(semi-harmonic spring, Eq. (6)) it exhibits force only if it is stretched. The vertical compo- 
nent of this force always pulls the cargo towards the MT. Thus, the mean height between 
the cargo and the MT (which determines how many motors can bind at maximum) results 
from the interplay between this force and thermal fluctuations of the bead. In contrast, if 
the motor also behaves like a harmonic spring when compressed, it once in a while may also 
push the cargo away from the MT. This results in less motors being close enough to the 
MT for binding than in the case of the cable-like behavior of the motor stalk. Consequently, 
given the same Zq and Ntot-, the cargo is on average transported longer distances when pulled 
by "cable-like motors". 

In order to apply the theoretical prediction for the mean run length of a spherical 
cargo particle, Eq. (19), we need to determine proper values for the three parameters 
flfe = nb/Ntot, U, Nm,max- From the simulations we measure the mean velocity of the sphere 
U . It turns out that U is up to 15 % less than the maximum motor velocity vq due to 
geometrical effects. Depending on the point where the motor is attached on the sphere, 
some motor steps may result mainly in a slight rotation of the sphere instead of transla- 
tional motion of the sphere's center of mass equal to the motor step length 5. For Nm^max 
we choose the overall measured maximum value from all simulation runs for given Ntot and 
polymer model of the motor. Then, we use the remaining parameter, the area fraction ab, as 
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Iq, motor-model 


fit value for at 


measured {h) — » ab{{h)) 


50nm, spring Eq. (5) 


0.00211 


7-14 nm 0.0039-0.0034 


50nm, cable Eq. (6) 


0.0026 


4-11 nm 0.006-0.0055 


80nm, spring Eq. (5) 


0.00403 


8-14 nm 0.0082-0.0076 


80nm, cable Eq. (6) 


0.00518 


4-11 nm 0.0085-0.0079 



TABLE II: Obtained fit values for the area fraction at for different Iq and the two applied polymer 
models. For comparison the area fraction which is obtained from the measured mean distance 
(h) is also displayed, (h) is measured for fixed Ntot, the left boundary of the provided interval 
corresponds to the largest Ntot- 

a fit parameter to the numerical results. The fits are done using an implementation of the 
Marquardt-Levenberg algorithm from the Numerical Recipes [46]. The resulting values 
are summarized in Tab. II. The theoretical curves for those parameter values are shown 
(dashed lines) in Fig. 7a. The increase in for larger resting length and the cable model is 
in excellent accordance with the above discussed expectation. An independent estimate of 
the area fraction can be obtained by measuring the mean distance (h) between cargo and 
MT and calculating the area fraction ab as = ab{{h)). For comparison those values are also 
given in Tab. II. They turn out to be about 60 % larger than the values for at obtained from 
the fit. This indicates that the height of the sphere above the MT (and therefore also the 
area fraction) is a fluctuating quantity that is not strongly peaked around some mean value. 
Then, because of the non-linear dependence of ab on h we have in general {ab{h)) ^ ab{{h)). 

Fig. 7b shows the mean number of bound motors (the average is obtained over all 
simulation runs) as a function of Ntot (symbols with error bars). The lines in Fig. 7b are 
plots of Eq. (20) using the same parameters as for the correspondig lines in Fig. 7a. One 
recognizes that again the theoretical prediction and the measured values match quite well. 
This means that on the level of the mean run length and mean number of bound motors 
the Poission average that was introduced in Sec. II F works quite well, even though the 
number of motors that are in range to the MT is not a constant during one simulation run 
(cf. Sec. IIIC). The large error bars for the cable model data in comparison to the spring 
model results partly from a poorer statistics (for the Iq = 80 nm cable simulations the 
number of runs is in the range of some hundreds only). In addition, for cable-like motors the 
fluctuations in the area fraction are much larger than for spring-like motors as repulsive 
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FIG. 8: Combination of data from Fig. 7a,b for the resting lengths Iq = 50,65,80 nm. (Axf,) is 
shown as a function of A'^;,. For the dashed hne (theory) Eq. (19) and Eq. (20) were combined 
with a truncation of the sums at Nm,max = 18. (Parameters: Had = 5 s, eo = 1 s, Ag = 125 s, 
At = 10-Veo, N ~ 10^.) 

spring forces stabilize the distance between the sphere and the MT. Therefore, the width of 
the distribution density of the number of bound motors is larger for the cable model than 
for the spring model. Fig. 7b also shows that for a bead radius of 1 /im, around 80 motors 
have to be attached to the bead, otherwise binding and motor stepping become unstable 
because there are less than two motors left in the binding range. 

In Fig. 8 the simulation results of Fig. 7a,b are combined into one plot. Here we show 
the measured mean run length (Ax;,) as a function of the mean number of bound motors 
Nb. All curves collapse on one master curve that can be parametrized by rtb = abNtot, he., 
the product of the fit parameter af, and the total number of motors on the sphere. The 
fact that all data points turn out to lie on one master curve again demonstrates the good 
applicability of the theoretical predictions to the simulation results. The curve shown in 
Fig. 8 has a positive curvature in the semi- logarithmic plot. This turns out to be an effect 
of the finite truncation of the sums in Eq. (19) and Eq. (20). 

C. Distribution of motors in binding range 

For the evaluation of the numerical data in Sec. IIIB we assumed that the number of 
motors that are in binding range to the MT are constant in time and that this number is 
drawn from a Poisson distribution for every individual run. We now further examine this 
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FIG. 9: Histograms for the number of motors that are in binding range to the MT. Symbols refer 
to simulation results for different values of the total number of motors on the sphere Ntot- Lines are 
Poisson distributions with mean value that is propotional to Ntot- (a) Resting length Iq = 50 nm, 
spring model, Eq. (5). (b) Resting length Iq = 50 nm, cable model, Eq. (6). (Parameters: -Kad = 
5, eo = 1, A° = 125, At = 10-Veo, N = 2- 10"^.) 

aspect in order to see to what extend this assumption is fulfilled. First, we measure directly 
the distribution of the number of motors umt in binding range to the MT. For this we 
count Umt at every numerical time step during one simulation run and repeat this for a 
large ensemble of runs (A^ ~ 10^). Thus, the histograms (i.e., approximately the probability 
distributions we are looking for) that are obtained in this way for the relative frequency of 
Umt are not based on the assumption of constant Umt for every single run. 

Fig. 9 shows examples of such histograms (symbols) for a series of different values of the 
total number of motor coverage Nt^t- For Fig. 9a we used the spring model for the motor 
polymer, Eq. (5), and for Fig. 9b the cable model, Eq. (6). In both cases the resting length 
of the motor protein is Iq = 50 nm. In addition. Fig. 9 also displays probability distributions 
(solid lines) that are obtained from the Poission distribution given that at least one motor 
is in binding range, i.e., p{n) = fi"'e~'^ /n\{l — e~"), with mean value /i given by /i = ^qNioi 
that can be parametrized by some variable /iq. The parameter /io was chosen to be 0.015 and 
0.0165 for the spring and cable model, respectively, by matching the Poisson distribution 
to the simulation data. For the spring model (Fig. 9a), the fit is excellent for all values of 
Ntot- One must note however that these distributions are not given by Eq. (18) as indicated 
by the fact that the parameter /io is much larger than the area fraction determined in 
the previous section. Instead one needs to account for the correct weighting factors from 
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FIG. 10: The relative frequencies of the number of motors that are in binding range to the MT 
during a single run are shown for six different sample runs. For the motor protein the cable model 
with resting length Iq = 50 nm is used. (Parameters: Ntot = 400, iTad = 5 s, eo = 1 s, = 125 s, 
At = 10-Veo). 

the different run lengths. If for the moment we consider the number hmt to be a constant 
during one run, then the distribution Eq. (21) can be considered to be a useful estimate for 
the data in Fig. 9. Taking the limit Nm,max — >■ oo in Eq. (21) we obtain 

(afc(l + W6o))"«^-a^^^" e-"^ ^ ^23) 

with h = af,(l + TTad/^o)- Thus, the result is approximately — except for very small timt = 
1, 2 — again a Poisson distribution with parameter b. With iTad/^o = 5 and at = 0.0021 from 
Tab. II we get b = 0.013 which is very close to the value fiQ = 0.015 used in Fig. 9a. 

For the cable model (Fig. 9b), the Poisson fit using a single value for /iq works well only 
for the smaller values of Nfot- For large Ntot the data points cannot be fitted by a Poisson 
distribution. It rather turns out that the ratio of mean value and standard deviation becomes 
less than one in these cases. 

Finally, we shall note that despite the good agreement between the simulation data and 
the the estimate from Eq. (23), which was based on the assumption that Umt is constant 
during one run, umt is in fact not constant, but a fluctuating quantity. The fluctuations 
result partly from thermal fluctuations of the height and orientation of the sphere and partly 
from orientation changes of the sphere that are induced by motor forces. Fig. 10 shows a few 
sample histograms for the frequency that umt motors are in binding range to the MT, i.e. 
either bound to the MT or unbound within the binding range, during a single run. These 
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FIG. 11: Measured probability distribution density for the escape rate e given e > eq. The data 
was obtained for different values of Nfot- For the motor proteins the full spring model, Eq. (5), 
with resing length Iq = 50 nm was used. (Parameters: TTad = 5, eo = 1, = 125, At = 10~^/eo, 
= 2 • 10^.) 

examples clearly show that hmt takes different values during one run. 
D. Escape rate distributions 

Diffusive motion of the cargo sphere directly influences the length of the pulling motor 
proteins and therefore also the escape rate e. The dependence of the escape rate on the 
motor length x for x > Iq is obtained by inserting the polymer model Eq. (6) and Eq. (5), 
respectively, into the Bell equation, Eq. (8). For x < Iq the escape rate is given by e = eg. 
For low viscous friction we assume x to be distributed according to a Gaussian distribution. 
Then, we expect the probability distribution density p(e) for the escape rate e to be given 
by a log-normal distribution density 

2kBTK* J' ^ ^ 

Here, e denotes the escape rate associated with the mean motor length in the extended state, 
K* is an effective spring constant that depends, e. g. on the number of pulling motors, and 
6 is a normalization constant that is obtained from the condition 

/■oo ^ 

/ p{e)de = 1. 
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exp 
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Fig. 11 shows the measured probabihty density for e > eo and different values of Ntot- It 
turns out that the log-normal distribution in Eq. (24) matches well to the measured data 
(not shown). Fitting Eq. (24) to the data for the effective spring constant k* it turns out 
that K* increases with increasing Ntot- This makes sense and illustrates that for several 
motors pulling the motors behave as a parallel cluster of springs. 

The agreement of Eq. (24) with the simulation data for the escape rate distributions 
suggests that the main source of force acting on the motors and increasing the unbinding 
rate is due to thermal fluctuations of the micron-sized bead. This is different from what has 
been reported in experiments with nano-scaled two-motor complexes [16]. There it has been 
argued that the forces between the two motors arising from their stochastic stepping lead 
to an increased unbinding rate of these motors. 

E. Cargo transport against high viscous friction 

Except for the single motor simulations of Sec. Ill A, all simulation data discussed so 
far were obtained for a viscosity of 1 mPa s, corresponding to a water-like solution. We 
mentioned in Sec. II B that load sharing between several motors may lead to cooperative 
effects at high viscous friction. We now analyze this further by performing simulations at 
viscosities much larger than that of water (i.e., when the viscous friction on the bead is 
comparable to the internal friction of the motor protein). To do this we need to measure the 
velocity of the bead depending on the number of pulling motors. Because of the nature of 
the stochastic process describing the position of the cargo the instantaneous cargo velocity 
is however not well defined [47] . Therefore, in order to measure the cargo velocity U we have 
to average over some time interval At. If no motors were pulling the velocity distribution 
density is given by a Gaussian, 



with diffusion constant D = ksTafi^r^^ (cf. Sec. II A). Furthermore we assume a constant 
height of the sphere so that the mobility coefficient /i"^ is a constant in time. The width of 
the distribution density, Eq. (25), is the smaller the larger At is. So in order to suppress 
fluctuation effects it seems appropriate to average over a large time interval At. On the other 
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FIG. 12: (a) Propability density of the cargo particle's velocity U that is obtained by averaging 
over a time interval of At = 0.02 s for different values of the viscosity 77. The inset shows the 
mean velocity as a function of the r]. (b) The mean velocity of the bead given that n motors 
are simultaneously pulling is plotted as a function of n (symbols). For comparison the theoretical 
expectation according to Eq. (13) is plotted, too (lines), (c) For high viscosity rj = 100 mPas the 
conditional velocity distribution density given that a certain number of motors is pulling is plotted. 
(The full harmonic spring model was used; parameters: = 80 nm, Ntot = 200, numerical time 
step At = 10~^, other parameters as in Tab. I.) 



hand the number of pulling motors changes with time because of binding and unbinding. 
Thus, in order to measure the velocity given a certain number of motors At should not be 
too large in order to get enough such events. Here we choose At = 0.02 s which corresponds 
to a typical camera resolution of 50 Hz. 

Fig. 12a shows the measured velocity distributions for three different values of the viscos- 
ity, ?7 = 1, 10, 100 mPa s. In the inset of Fig. 12a the mean velocity is plotted as a function of 
the viscosity. It turns out that shifting the distribution Eq. (25) by the corresponding mean 
velocity, the single peaked function p{U, At) fits quahtatively well to the distribution shown 
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in Fig. 12a, especially the dependence of the width of the distributions on r] is correctly 
predicted by Eq. (25). Also a decrease of the mean velocity of the cargo particle is observed 
with increasing viscosity resulting from the increased frictional load. It must be emphasized 
that only a single peak is observed in Fig. 12a, even though the bead velocity is expected 
to depend on the number of pulling motors, which should lead to multiple distinct peaks 
[7] . One reason that we do not observe multiple peaks is the small value of the time interval 
At, which leads to broad peaks with a peak width governed by diffusion of the bead (cf. 
Eq. (25)) and thus makes it impossible to separate different peaks. Using larger values of 
At leads to smaller peak widths, but also to poorer statistics as less measurement points are 
obtained, so that again distinct peaks cannot be resolved. Even if we do not use a constant 
time interval, but average over the variable time intervals between two subsequent changes 
in the number of bound motors [9], distinct peaks are very hard to separate (not shown). 
This does however not mean that the bead velocity is independent of the number of pulling 
motors. Indeed if we plot the conditional velocity distribution calculated over all intervals 
in which the bead is pulled by a certain fixed number of motors, we see a clear shift in the 
average velocity (Fig. 12c). This shift is however masked by the width of the distributions 
in Fig. 12a. 

In Fig. 12b the average of all measured velocities given that exactly n motors are pulling 
is plotted as a function of n, again for the three different viscosities r] = 1, 10, 100 mPa s. 
For rj = 1 mPa s the viscous friction for the bead is about l/f4!:x ~ ^ ■ 10"^ Ns/m. The 
internal friction of the motor is 1/fim = Fs/vq ^ 6 ■ 10"^ Ns/m, i.e., about two orders of 
magnitude larger than 1/ fi^J:^- According to the analysis at the end of Sec. II B we therefore 
expect that the mean velocity is independent of n if all motors equally share the load. The 
numerical data in Fig. 12b shows that the mean velocity exhibits a weak dependence on n 
with a maximum at about n = 4,5. At the higher viscosities the numerical results show 
that the mean bead velocity increases with increasing n, which indicates that the motors 
share the load. The simulation data however deviate clearly from the estimate given by 
Eq. (13), which is indicated by the lines in Fig. 12b. This discrepancy indicates that the 
load is not shared equally among the motors or that only a subset of the bound motors are 
actually pulling the bead. Besides geometrical effects one reason why this is the case is that 
the escape rate eo is rather high and at high frictional load is even further increased making 
the lifetimes of motors in the pull state rather short. Then, if a new motor binds to the MT 
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FIG. 13: Frequencies of the motor forces in walking direction relative to the total load force on 
the cargo particle for different numbers of pulling motors, (a) Viscosity rj = 1 mPa s, escape rate 
eo = 1 s-^ (b) t] = 100 mPa s, eo = 1 s~^ (c) t] = 100 mPa s, eo = 0.1 s~^ (d) rj = 1000 mPa s, 
eo = 0.1 s^^. For the other parameters the same values as in Fig. 12 were used. 



often another motor detaches already before a stationary state is reached in which all the 
motors equally share the load. 

In order to investigate the last point in more detail we consider explicitely how the force 
is typically distributed among the pulling motors. For this we count the number n of motors 
attached at each numerical time step and measure the force experienced by each motor in 
the direction along the microtubule. For a given number n, n such motor forces can be 
measured, Fm]x,i = 1, . . . ,n. To suppress effects from fluctuations in the overall load we 
then calculate the reduced forces /„ := Fm]x/ Y^'^=i-^in,x- Given the histograms for these 
quantities measured over many time steps and simulation runs we obtain an approximation 
for the probability distribution density of the relative load of the motors. Fig. 13 shows 
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results for such histograms obtained at different viscosities and two different values of the 
unstressed escape rate eo- Fig. 13a shows the results for 77 = 1 mPa s and eo = 1 s~^. 
One can see some symmetries that result from the definition of the reduced forces. This 
is especially emphasized for the case of n = 2 pulling motors. The distribution density of 
the reduced forces has a mirror symmetry about the value /2 = 0.5. Besides this artefact it 
turns out that for n > 2 the distribution densities are strongly peaked at zero force. Thus, 
for low viscous friction the overall load results mainly from fluctuations. Such loads are 
typically experienced by a single motor only, whereas the remaining motors are more or less 
un-stretched. In Fig. 13b the viscosity is 100 times larger than that of water which causes 
an appreciable load on the pulling motors. However, load sharing effects are hardly visible 
except for the maxima around /2 = 1/2 and /s = 1/3, which are rather broad and thus 
hard to resolve. On the other hand there is still a strongly peaked maximum at /„ = 0. 
This somewhat surprising observation is due to the rather high escape rate, which is even 
increased by the load force (cf. Eq. (8)). Therefore, the binding time of the motors is 
shortend. On the other hand motors that bind to the MT are initially unstressed (i. e., carry 
zero load) in our model. Thus they always contribute to the /« = peak and may already 
escape from the MT before the load is shared equally by all pulling motors. When the 
escape rate is reduced to eo = 0.1 s"^ as done for Fig. 13c, d clearly visible maxima around 
/„ = 1/n appear in the histograms that indicate that the load is equally shared by the active 
motors. In Fig. 13d where we used the extremely high value of 10^ mPa s for the viscosity, 
these peaks are very pronounced. The arrows in Fig. 13c, d indicate the relative force values 
1/z, z = 2, 3, . . .. It turns out that the peaks are not exactly located at these values which is 
again due to the binding and unbinding process of the motors. 

In summary, we found that at high loads the pulling motors tend to arrange in such a 
way that the total load is equally shared amongst them. However, for typical escape rate 
values of kinesin-like motors this process often takes more time than the lifetime of a state 
of a certain number of motors bound to the MT lasts, thus preventing cooperativity in the 
sense of equal load sharing. 
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IV. DISCUSSION AND OUTLOOK 



The main purpose of this paper is to introduce a novel algorithm called adhesive motor 
dynamics as a means to study the details of motor-mediated cargo transport. Our algo- 
rithm is an extension of existing adhesive dynamics algorithms developed to understand the 
physics of rolling adhesion [22, 26]. Basically our method allows to simulate the motion 
of a sphere above a wall including hydrodynamic interactions and diffusive motion by nu- 
merically integrating the Langevin equation, Eq. (3). In addition, motor-specific reactions 
such as binding to the microtubule and stepping are modeled as Poisson processes and then 
translated into algorithmic rules. The parameters and properties by which the motors are 
modeled are based on results of single-molecule experiments with conventional kinesin. A 
first favorable test for the algorithm was provided by measuring the force velocity-relation 
at different viscosities and external load forces and by comparing the results to the input 
data as done in Sec. Ill A. 

Next we measured the run length and the mean number of bound motors as a function 
of the total number Ntot of motors attached to the sphere. The same quantities have been 
previously calculated based on a one-step master equation model [7]. However, this has 
been done as a function of a fixed number of motors that are in binding range to the 
microtubule. In practise and also in our simulations, this quantity varies in time. In Sec. II F 
we Poisson-averaged the theoretical predictions thus rendering it possible to compare theory 
and simulation results. Using the area fraction on the cargo from which the microtubule is 
in binding range for the motors as a fit parameter, we found good agreement between theory 
and simulations for both the mean run length and the mean number of bound motors. Note 
that the latter one cannot be measured in typical bead assay experiments. 

We also determined the mean separation height between cargo and microtubule and 
found (h) = 4 — 14 nm. Modeling the motor stalk as a cable resulted in smaller distances 
than using a full spring model for the motor stalk. A recent experimental study using 
fluorescence interference contrast microscopy found that kinesin holds its cargo about 17 nm 
away from the MT [37]. Our smaller distance probably results from neglecting any kind 
of volume extension (except binding site occupation) of the motor protein, the simplified 
force extension relation applied to model the stalk behavior, and neglecting electrostatic 
repulsions. These effects, however, could easily be included into our algorithm, e. g., by using 



33 



hard core interactions that account for the finite volume of the motor protein segments. 

In Sec. Ill C we explicitely demonstrated that the theoretical assumption of having a fixed 
number of motors in binding range during one walk is not justified (cf. Fig. 10). Nevertheless 
the theoretical results agree well with the simulations. This might be explainable by the 
observation that averaged over many runs the distribution of motors in binding range appears 
to be Poissonian (cf. Fig. 9). Thus, on the one hand fast fiuctuations in the number of motors 
in binding range around some mean value are not visible. On the other hand periods in which 
this number fiuctuates around the same mean value can be treated as a complete run. Thus, 
averaging over these smaller runs (i.e., which end after the sphere was e.g. rotated visibly 
and not after the last motor unbinds) has the same effect as averaging over complete runs 
(i. e., which end after the last motor unbinds). 

An interesting question is to what extend several motors can cooperate by sharing load. 
We have addressed this question for the case of several motors pulling a cargo particle 
against high viscous friction. One of the advantages of our algorithm is that we can measure 
the velocity of the bead and at the same time also the number of simultaneously pulling 
motors. Thus, in Sec. Ill E we tried to check whether the explanation of Ref. [9] is correct also 
under the assumptions of our model, especially for the parameter values given in Tab. I. Our 
simulations show that the speed of the cargo increases with the number of pulling motors for 
high viscous friction, in agreement with experimental results [48]. Our simulations however 
show pronounced deviations from the quantitative predictions based on the assumption 
that the load is shared equally among the bound motors. Furthermore, as the average 
life time of a state with a certain number of pulling motors is rather short the different 
velocities expected for different numbers of instantaneously pulling motors were smeared 
out by diffusion. Similarly when directly measuring how the total load force is distributed 
to the different motors pulling, no equal load sharing could be observed for the escape rate 
of about 1 s~^. We observed equal load sharing only when we used a ten- fold smaller escape 
rate, in order to increase the life time of the motors in the bound state. 

Another interesting question in this context is whether the velocity distribution exhibits 
several maxima if the cargo is pulled against a viscous load, as observed in several experi- 
ments in vivo [9, 49]. For example. Hill et al. [9] found that vesicle in neurites move with 
constant velocity for some period of time and then switch to another constant velocity in a 
step-like fashion. The distribution of velocities (measured over time intervals of the order of 
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1 s) was found to have peaked at integer multiples of the minimal observed velocity. These 
peaks were interpreted as corresponding to different numbers of simultaneously pulling mo- 
tors, which equally share the visoelastic load excerted by the cytoplasm [9] (cf. also Eq. (13)). 
Indeed, both an earlier model for motor cooperation [7] and our present description predict 
that equal sharing of a large viscous load leads to such a velocity distribution. In our simu- 
lations, we could however not resolve multiple peaks, presumably because the peaks are too 
broad to be resolved. The latter results from a combination of the way how we measure the 
velocity and from the fast dynamics of motor unbinding as discussed in Sec. Ill E. 

As already mentioned above, the framework of our method is rather general. Therefore 
various model variations can be easily implemented and probed in simulations. Here, for 
example we modeled the motor stalk by two versions of a simple harmonic spring: the 
cable model, which represents a molecule with an intrinsic hinge, and the spring model. 
More advanced force-extension relations could easily be incorporate in Eq. (4) in order to 
probe the influence of more realistic polymer models on the transport process. Similarly, 
the force dependence in unbinding from the microtubule and stepping can be altered to 
explore the impact onto macroscopic observables like the mean run length or the speed 
of the cargo. Furthermore, the algorithm could easily be adapted to study beads to which 
clusters of motors or defined motor complexes (such as those in ref. [16]) are attached. Thus, 
the algorithm described in this paper can be regarded as a link between purely theoretical 
models and data from in vitro experiments. 

Another interesting question for future applications of our method is how cargo transport 
works against some external shear flow. Since our model is based on a hydrodynamic 
description, flow can easily be included in the dynamics of our model. For these studies 
the Langevin-equation, Eq. (1), has to be extended by additional terms accounting for the 
effect of an incident shear flow [23, 28]. Then, two opposing effects exist characterized by 
the step rate and the strength of the shear flow, respectively. Their interplay together with 
the rates for binding and unbinding Tiad and e, respectively, determine whether the cargo 
moves in walking direction or in flow direction. Experimentally, such a setup might provide 
interesting perspectives for biomimetic transport in microfiuidic devices. 
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APPENDIX A: ADHESIVE MOTOR DYNAMICS 

The Langevin equation, Eq. (3), and the motor dynamics rules explained in Sec. II B are 
connected by the following algorithmic rules that apply in each time step At: 

(i) The sphere's position and orientation is updated according to Eq. (3) (for an explicit 
description see Ref. [26]). 

(ii) The positions where the motors are attached to the sphere in the flow chamber coor- 
dinate system are calculated. 

(iii) For each motor that is bound to the MT its load force is calculated. Then stepping 
is checked according to the stepping rate derived from Eq. (10). If the motor steps 
forward the load force is again calculated as motor length/direction has changed. 

(iv) For each motor that is not bound to the MT binding is checked according to the 
procedure explained in the main text (Sec. II B). 

(v) For each active motor (i. e., bound to the MT), the contribution to is calculated. 

(vi) Each motor that is bound to the MT can unbind with escape rate e given by the Bell 
equation, Eq. (8). 

A motor that escaped from the MT in one time step can rebind to the MT in the next time 
step according to rule (iv). The same Monte-Carlo technique that is explained in the main 
text (Sec. II B) to decide whether binding occurs or not is also used for the decission on step- 
ping and unbinding. For the pseudo random number generator we used an implementation 
of the Mersenne Twister [50]. 
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APPENDIX B: MOTOR DISTRIBUTION ALGORITHM 



Initially the center of the sphere is located at position (0, 0, R + lo + hMr) directly above a 
microtubule binding site (cf. Fig. 1). The first motor that is distributed is initially fixed at 
position (0, 0, /q + ^a/t) (relative to the chamber coordinate system) with its tail. The head 
is bound to the microtubule at (0, 0, Hmt)- Thus the initial distance between the motor and 
the microtubule is given by the motor resting length Iq. For the distribution of the other 
Ntot — 1 motors on the sphere we use an hard disk overlap algorithm similarly to the one 
that was described in Ref. [22]. For each of these motors two random variables are chosen 
ri from the uniform distribution on (0, 27r) and r2 from the uniform distribution on (0, 1), 
respectively. Then, the motor is located on the sphere's surface at the spherical coordinates 
(ri,arccos(l — 2r2)) and possible overlap to already distributed motors is checked. If no 
other motor is located within a ball of radius O.I/q around the just distributed motor, then 
its position is kept, otherwise a new pair of random coordinates are drawn until no overlap 
with other motors exists. 
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